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2 Introducción 


Resolver una ecuación de tipo f(x) = 0 no siempre es fácil, es más, si se tratara de un sistema 
de ecuaciones de la forma: 


fi (Union) =( 
fa yt) = 0 
de (E 0 


, el problema se complica mucho más. Sin embargo hay métodos numéricos que posibilitan 
encontrar una aproximación aceptable de la solución de dichos sistemas, en este material se 
busca exponer métodos numéricos que logren conseguir soluciones para sistemas de ecua- 
ciones no lineales (SENL), en algunos casos estos métodos son análogos a los utilizados en 
el caso lineal como lo son los métodos de puto fijo y de Newton, además de un método 
cuasi-Newton como lo es el método de Broyden que propone una mejora en el método de 
Newton relajando algunos cálculos del proceso en el algoritmo de Newton, se expone también 
del método de descenso más rápido que aunque no es tan exacto, puede dar resultados que 
sirvan como aproximación inicial para otros métodos y para finalizar, se introduce el método 
de continuación que utiliza un parámetro auxiliar para alcanzar el valor deseado de la aprox- 
imación. 


Para cada uno de estos métodos se propone un código correspondiente en lenguaje Oc- 
TAVE/MATLAB y al final de este material se resuelve un problema de aplicación por medio 
de la implementación de uno de estos métodos, se omiten las demostraciones de los teoremas 
enunciados pero se hace énfasis en la motivación del algoritmo asociado a cada método. 


3 Preliminares 


Un sistema de ecuaciones no lineales con n variables y n ecuaciones tiene esta forma: 


los) =0 
ToDo 1 L) =0 
sa) =0 
: donde cada función f;, para cada ¿=1,2,...,n está definida de esta manera: 
FR” — R 
(11,L0,.., En) > filo1,%a,..., En) 


: donde alguna función es no lineal respecto de alguna variable x,, cada función f; se puede 
llamar campo escalar o función escalar. 


Luego se puede representar el sistema, definiendo una función F de esta manera: 


F:R" —> R” 


AS 0 TT Cf A A A A 


Bajo ciertas condiciones esta función se llama campo vectorial, pero para propósitos del tema 
se hace énfasis en el uso de negrita para denotar este tipo de funciones, así a menos que 
se indique lo contrario, estas notaciones por ejemplo: F, G y H denotan funciones con las 
características anteriores. 


Te 7 E 7 t 

Para facilitar cálculos se introduce la notación x = (11,%3,...,Ty) para denotar el vector 
cuyas componentes son 1,,%2,...,T,, notar que se usa negrita, así por ejemplo y,s y p son 
vectores en |R”. 


Para medir” la distancia entre elementos (vectores) de R” se introduce lo siguiente: 


e Definición 1: La norma infinito de x = (%1,%2,...,tn) se denota |lx|l.o. y se define así: 
ix], = max [5] 
1<i< 


Luego la distancia entre dos vectores, x y y por ejemplo se define como la norma infinito de su 
diferencia así: ||[x — y|lo, , esto será muy útil en el momento de comparar las aproximaciones 
que se vayan generando. 


4 Métodos numéricos para la solución de sistemas de 
ecuaciones no lineales 


4.1 Método de punto fijo 


e Definición 2: Una función G: D CR” —> KR” tiene un punto fijo: p € D si G(D) = D 


e Teorema 3: Sea D = laj,b1] x [az,ba] x --- x [a,,b,] para conjunto de constantes 
41,092,...,0n y b1,b2,...,by. Sea G: D CIR” —> R” un función continua con la propiedad 
de que G(x) € D siempre que x € D (G tiene un punto fijo en D). Supongamos que las 
funciones componentes de G tienen todas sus derivadas parciales continuas y que existe una 
constante K < 1 tal que: 


'9) ¡AX K e 
| gx) < — siempre que x€ D 
Ox; n 
: esto es para toda ¿ =1,2,...,n y toda función componente g, con ¿=1,2,...,n. 


Entonces la sucesión ENYA definida desde alguna x% € D seleccionada arbitrariamente 
y generada por x*) = G (xD), para cada k > 1 converge al punto fijo: p € D de G y es 
único, además: 


k 
[px po, < E [pa 00) 


[0,0] 


4.1.1 Algoritmo 


xd eDCcR” 
x= G (xD) ,%> 1 


: donde G cumple las condiciones del teorema 3. 


4.1.2 Ejemplo 


Demostrar que el sistema : 
xíi—101+12+8  =0 
xix + x1—10x2+8 =0 


, tiene solución única en D = [a a) [O < 21,19 € 1.5) y se calcular dicha solución con 
Y=:(0,0 


? 
con aproximación inicial x( )' y tolerancia de 107?. 


o Solución: 


Se plantea el problema como una búsqueda de punto fijo, definiendo una función G : R? —> 
R? donde sus funciones componentes q, : R? — R y g2 : R? — R se obtienen despejando 
las ecuaciones del sistema así: 


2 2 
za  X Da O 
91 (21, 22) = Li 7 
2 
8 
(11,19) = 22= A 


Notar que G(D) € D porque 5 € ca < 1.25 y también $ < pp < 1.2875 


Obteniendo las derivadas A de G, se tiene que: 


Oglz)| _ |9g()] _, 
Óx1 0x2 
Ogi(x) _ 22 < 1.5 03 Oga(x) a! Z 15 _ 0.3 
Ox, 5 9 0x1 9 y 
09, e 
>| aa) <0.3 ;paratoda 1,j=1,2,...,n y toda € D 
J 


Haciendo E = 0.3, se tiene que K = 0.6 < 1, se nota además que cada derivada parcial es 
continua en D y aplicando el teorema se tiene que G posee un único punto fijo en D, luego 
aplicando el algoritmo para x(% = (0, 0), se tiene para k > 1: x*) = (21,24) = G(x*-D), 


(k) (%) 


donde los valores de 1; ” y 13” vienen dados así: 
2 2 

o A) + 

1 
10 
Ro) (ny? (1 

a Y 

Lo = 


> Primera iteración: 

x0) = (0,0) > xU =G (x0) = (0.8, 0.8)' 

> Segunda iteración: 

x(Y = (0.8,0.8)* > x%) =G (xW) = (0.928, 0.9312)! 
> Iteraciones siguientes: 


Tabla 1: Iteraciones de método de punto fijo 


k 41) 5) [e — (0 ¡0 
1 Cuco ) (ome ) [0800000 
2 (oso ) (ossrzmo ) | 0191200 
, omo a 0.0448317 
4 ios! ta 0.0165339 
TEME 
6 | Cogosrasr ) (omo) 00025362 


MEME 
Casos) Comme ) | osos 
o (usaras) Coiimasas ) | omoron 
ol (eomios ) | (ommase) | oso 
dl al mica 0.0000258 
12 oa Ea 0.0000103 
| (esmas) (esemar) | 0oomos 


Se considera la aproximación: x("%) = (0.9999973, 0.9999973)' con tolerancia de 107?, el valor 


real del punto fijo es x = (1, 1). 


4.1.3 Código 
Código 1: Implementación del método de punto fijo 


function [x i]= PuntoFijo_SENL(Fsis, x0, TOL, n=35) 
AINPUT: 
/Fsis: función de IR” a R” asociada al sistema 
40: aprorimación inicial 
/TOL: tolerancia permitida 
/n: numero máximo de iteraciones (35 por defecto) 
AOUTPUT: 
4x3: aprorimación final 
/4t: indice de la ultima iteración 
fprintf£('Wtt>> Método de Punto Fijo'); 
fprintt (INtx (NN O NN RRA) A oo Nao 
o A A A A 
O A 

x = Fsis(x0); 


err = max(abs(x-x0));%para comparar con la tolerancia (es la norma infinito) 


fprintf ('pdXt/10.7£Xt410.71Xt/10.7fWn',1,x0(1),x(1),err); 
for j= 2:length(x) 
fprinti Ny 0 EN AO EN CA 0) a 
endfor 
ca A A A A A O A A Na, 


if err < TOL /se comprueba sí ya se satisface la tolerancia 
return 


else 


x0 = x; / sí no se satisface, se realiza el mismo proceso 
Arefrescando el valor de x0 
endif 
endfor 
/se llega al máximo de iteraciones 
disp("Numero máximo de iteraciones alcanzado") 
return 


endfunction 


4.2 Método de Newton 


El algoritmo del método de Newton para una sola ecuación no lineal viene dado por: 


IN EA CA, 


PD) 


Para asociar este resultado a funciones de |R” se se introduce la siguiente definición: 


e Definición 4: Sea F : R” > ¡R” una función tal que sus derivadas parciales de primer 
orden existen en RR”. La matriz Jacobiana de F denotada por Jp está definida como una 
matriz de tamaño n Xx n que tiene esta forma: 


f Of Of 
Ofa 942 (x 2B. (y 
Je (x) — dd de de 
(9) n (2) n (2) n 
92 (x) 9 (x) e 32 (x) 


Ahora en R” se puede introducir la siguiente fórmula para el cálculo de las aproximaciones: 
x 0) = xD — [jp (x4-0)] > F (241) 


, donde la ecuación matricial anterior tendrá solución siempre que el det (Jp (xD) sea no 


nulo para cada k > 0, notar que para cada iteración es necesario calcular (Jp (x0=D)] DE 


4.2.1 Algoritmo 


x(0) e R” 
x() = xD) — [Jp (x4-D)]7 PE k=>1 
: donde Jp es la matriz Jacobiana de F y la ecuación tendrá solución siempre que det (Jp RA 
sea no nulo. 
4.2.2 Ejemplo 
Resolver el siguiente sistema con aproximación inicial x% = (0, 0,0)? y tolerancia de 107*: 


61, — 2cos (12123) — 1 =0) 
922 + y 1? +sen (13) +1.06+0.9 =0 
60x3 + 3e “12 + 107 — 3 =0 
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o Solución: 


F(x) = (6 — 2cos (12123) — 1, 912 + ya + sen (13) + 1.06 + 0.9, 
6013 + 32% + 107 — 3) 
6 2x3sen (1213) 2x2sen (1213) 
aj 21 9 0.5cos(13) 
r(x) E V13+sen(m3)+1.06 x?+sen(23)+1.06 
—3e UU go —3e "142 g 60 
> Primera iteración: 
F (x(%) = (3.0000, 1.9296, 31.4159)' 
6.0 0 as 0.1667 0 0 
Jp(x%) =| 0 9 0.4856 | > [Jp (x0)] *= O 0.1111 —0.0009 
0 0 60 0 0 0.0167 


=> x) = x0 — [fp (x0)] 7 F (x0) = (0.5000000, —0.1861424, —0.5235988)' 
> Segunda iteración: 


F (x)) =  (0.0094917,0.12472, 0.29262)' 


6 —0.1019 0.036227 
Jp (x) = | 0.55556 9 0.48113 
0.6129 —1.6463 60 
, 0.16648  0.0019006 0.000085278 
> [Jp (xU)] *= | —0.010171  0.11083  —0.00089488 


—0.0019797 0.0030217  0.016641 
=> x20 =x0 — [Jp (x0)] "E (x0) = (0.3306729, —0.2434028, 0.5194079)' 


> Iteraciones siguientes: 


Tabla 2: Iteraciones de método de Newton 


% TI) mo) [E 
0.0000000 0.5000000 

1 0.0000000 —0.1861424 0.5235988 
0.0000000 —0.5235988 
03333333 0.4981578 

2 —0.2426495 —0.1996068 0.0134645 
0.5235988 —0.5288264 
0.3306729 0.4981447 

3 —0.2434028 —0.1996059 0.0000131 
0.5194079 —0.5288260 
0.4981447 0.4981447 

4 | | —0.1996059 —0.1996059 0.0000000 

—0.5288260 —0.5288260 


Se considera la aproximación: xY = (0.4981447, —0.1996059, —0.5288260)' con tolerancia de 
10%, 


4.2.3 Código 


Código 2: Implementación del método de Newton 


function [x i]= Newton_SENL(Fsis, x0, TOL, n=35) 
AINPUT: 
/Fsis: función de R” a R” asociada al sistema 
4x0: aprorimación inicial 
ATOL: tolerancia permitida 
/n: numero máximo de iteraciones (35 por defecto) 
AOUTPUT: 
A/z: aprorimación final 
ht: índice de la ultima ¿iteración 
fprintf('Mtt>> Método de NewtonWn'); 
fprant ita ON EN ON lt) SON 
print AAA E A A A Nau 
Ol 
JO = jacob(Fsis, x0, TOL);%se calcula la matriz jacobiana aplicada a z0 
if det(JO) == O /se comprueba que sea invertible 
disp("Utilice otro método") 
return 
endif 
JI = inv(J0);%se calcula la inversa 
Fx = Fsis(x0);%se calcula el valor de F sobre z0 


x= x0 - JIxFx;%se aplica la fórmula de la aprorimación de Newton 
err = max(abs(x-x0));%para comparar con la tolerancia (es la norma infinito) 


%se refresca la tabla 
printt (ant OTE CONO 0 do ac eze 
for j = 2:length(x) 
fprint (UNION EN O AR 
endfor 
paar E ARA EE A A NoE 


if err < TOL %/se comprueba si ya se satisface la tolerancia 
return 
else 
x0 = x; / sí no se satisface, se realiza el mismo proceso 
Arefrescando el valor de x0 
endif 
endfor 
/se llega al máximo de iteraciones 
disp("Número máximo de iteraciones alcanzado") 


return 
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42 


endfunction 


Código 3: Aproximación de Jp en x 


function J = jacob(Fsis,x0,h) 

AINPUT: 

/Fsis: función de IR” a R” asociada al sistema 
40: aprorimación inicial 

4h: tamaño de paso 

AOUTPUT: 


4J: matriz jacobiana e F en el punto 20 


n = length(x0); dimensión 
J = zeros(n,n);%/se inicializa la matriz a llenar 


y = x0;%para las aprorimaciones de las derivadas 
Fx = Fsis(x0);%el valor de F en 20 
for i = 1:n/recorriendo filas 

fx = Fx(i);%se toma el valor requerido 

for j = 1:n/recorriendo columnas 


y(j) = y(G)+h;%se aumenta el correspondiente 


Fy = Fsis(y);%se calcula el valor de F en 20_j+h 
fy = Fy(i);%se toma el valor requerido 


J(1,j) = (fy-£x)/h;%aprozimación de la derivadas 
va x0;%para comenzar otra vez 
endfor 
endfor 


endfunction 


4.3 Método de Broyden 
Del método de Newton se tiene: x%) = x(4=1) — (Jp (x4-=D)] Hp (x=), 


El método de Broyden aplica el mismo algoritmo lo que lo convierte en un método cuasi- 
Newton que se diferencia del método de Newton al proponer una forma de aproximación para 
las matrices (Tp (x0-=D)] HS para cada iteración. 


Considerando f(x) = 0, un sistema de ecuaciones lineales entonces es de la forma Ax—b =0. 
Por el método de la secante se sabe que: 


6d (2) == (e 4=D) = a (20 — 0) 


Pero en el caso no lineal la igualdad anterior no se cumple aunque se puede lograr, eligiendo 
adecuadamente Az, tal que: 


F (200) —F (ac) — Ay (200) — el=1)) 


Se añade además la condición A¿(z) = Ax-1(z) siempre que (x0) — xD) z = (, es decir 


que para los vectores z ortogonales a (x0) — x 0-1) se tiene que Az(z) = Ax-1(Z) y con esto 
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cada Aj queda bien definida y se deduce de las condiciones anteriores: 


F (x4) -F (A 


di (50) — x (=D) 
F (x) -F (a 

= Ar-1 (4) — x (1) == Ar-1 

En F (xM) ==, F (xD) ==, Ar-1 (uedn == o 

SE as (x0) — x (6D) 

e (E (0) — FE (xD) — Az 1 (0 — x (4D) O o)! 

¡e (0%) — g (01) (x0) — g(—1)) 
Es decir: 


(Y — Ar-18r) ; 
2 k 
lIsrll 


: donde y, = F (xM) —F (0 y Sy =x 4 — x(4-D) 


Ar = Api A 


Luego se aplica el método de Newton deduciendo en cada iteración la matriz Az en sustitución 
de Jp y calculando su inversa. Para mejorar el método se plantea una forma de generar las 
matrices inversas directamente. 


e Teorema 5 (Fórmula de Sherman-Morrison): Si A es una matriz no singular y si x 
y y son vectores tales que y'4Ax 4 —1. Entonces A + xy? es no singular y además: 


E ATlxyiA7! 
A E =L 3% AT! NM A 
( + xy ) 1 + y Arlzx 
Aplicando la fórmula de Sherman-Morrison para A = Ap-1 , X = NS y y = Sk, Se 


tiene: : Ca 
(sy A Aj yr) S.A1 


EAT 
SAL 1Y k 


4.3.1 Algoritmo 


x0) e R” 
x0) = xD — ADE (xD) E>1 
(sy >= Aj 11) SAR 
SA 1Y 1 


ye =F (20) -E(x6-D) y 5, =x0 — x (00 


; donde Aj = [Je (x0)] 7, Az! = Ag, y 


4.3.2 Ejemplo 


Resolver el siguiente sistema con aproximación inicial x(% = (2, 2) 


In (12 + x2) — sen (2,12) =In(2)+In(r) 
ena 4 cos (1122) =.(0) 
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o Solución: 


F(x) = (In (2 + 23) — sen (2,23) — In(2) — In(m), 
en? + cos (1122)) 
2x1 A 2x2 ZE 
Ade 27 cos(1,13)12 27 cos(1113)11 
eri — sen (1119) 19 —e 172 — sen (1,12) 21 


> Primera iteración: 


F(x%) = (0.9984,0.3464)' 
1.8073 1.8073 0.1421 0.5 
= (0) = == 
AS Y [ets al % na El 


sy = Ay F (x0) = (-0.031211, —0.521105)' 
> x) =x0 — A7*F (x0) = (1.9686887, 1.4788947)' 


> Segunda iteración: 


F(xM) = (0.2638, 0.6583)' 


y =F (x) —F (x0%) = (-1.2621, -0.3120)' 


s; =xUY — x(% = (-0.031211, -0.521105)' 


0.1064 0.3300 


Aj' = Ay" + (1/84) 'y1) [(81 — 40"y1)5:40'] =| 03480 0.2623 


> x0) =xU — AUF (x%) = (1.7795004, 1.7433937)' 


> Iteraciones siguientes: 


Tabla 3: Iteraciones de método de Broyden 


> TI) mo) =D 
E onaooo lero, ) | A 
2 (rárososr ) | l rrasaos7 ) | 02004080 
(rasos ) | rimmsoosr ) | oos0ario 
rolas] | 004 
51 (armasesa ) | arrasa ) | 00000821 
o (rasa ) | rrasso9 ) | 00000303 
1 (rrmaass ) | irrassso ) | 00000010 
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Se considera la aproximación: x(") = (1.7724539, 1.7724539)* con tolerancia de 107?. 


4.3.3 Código 
Código 4: Implementación del método de Broyden 


function [x i]= Broyden_SENL(Fsis, x0, TOL, n=35) 

AINPUT: 

/Fsis: función de ¡R” a R” asociada al sistema 

40: aprorimación inicial 

ATOL: tolerancia permitida 

/n: numero maximo de iteraciones (35 por defecto) 
AOUTPUT: 

4x3: aprorimación final 

At: índice de la ultima iteración 

fprint£('Wtt>> Método de BroydenWn'); 

fprintf ('kWtx" (%-1)W6Xtx" (4) Wet | | x" (0) -x7 («-1) | [ooYn'); 
print IA IE A AE AAN A PO NoE 


JO = jacob(Fsis, x0, TOL);%se calcula la matriz jacobiana aplicada a TO 


v = Fsis(x0); 

if det(JO) == O /se comprueba que sea invertible 
disp("Utilice otro método") 
return 


endif 


A = inv(J0);%se calcula la inversa 
= (Arv).* (1); 
OS 


err = max(abs(x-x0)); 


if err < TOL %/se comprueba sí ya se satisface la tolerancia 


return 
else 
1 a 


priatiddNcA LONE OEA taa 0 a ferry 
for j = 2:length(x) 
fprintt (UNOS EN GALO TENEN ODER 
endfor 
print EA Wn') 
x0 = x;% si no se satisface, se continúa el mismo proceso 


endif 


for i = 2:n 


W = vv; 
v = Fsis(x0); 
y = vw; 
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44 


45 


46 


47 


48 


49 


50 


51 


52 


53 


54 


55 


56 


57 


58 


59 


60 


61 


62 


63 


64 


65 


66 


67 


68 


69 


70 


(Axy) .* (1); 


p= ((s')*z).x*(1); 

u= ((s')x*A); 

A = A+ ((s+z)xw).*(1/p); 
s = (Axv).*(-1); 

A OE SE 


err = max(abs(x-x0));%para comparar con la tolerancia (es la norma infinito) 


%se refresca la tabla 

print (Maty LO MEN ELO ENE LOS EN Odd erre 

for j = 2:length(x) 
fprintidNo 0 ENG IOEN AO 

endfor 

fprintíf(' 


if err < TOL %se comprueba sí ya se satisface la tolerancia 
return 
else 
x0 = x;% si no se satisface, se realiza el mismo proceso 
endif 
endfor 
/se llega al máximo de iteraciones 
disp("Número máximo de iteraciones alcanzado") 
return 


endfunction 


4.4 Método del descenso más rápido 


Este método determina un mínimo local para una función de varias variables de la forma 
g :|¡R” —>R. La relación entre un sistema de ecuaciones no lineales se debe a que un sistema 
lineal de la forma: 


fi (Lis Dos ari a) =0 
to (11, Ta, Aroa 6) =0 
A A 
, tiene solución x = (11, 79, ..., Tn)" justo cuando la función y definida como g(21,T3,..., Un) = 
(filas, ta, 2) + (falar, ta... En) +0 +(fn(t1, 22, ..., Ln)” tiene el valor mínimo. 
El algoritmo se describe de forma sencilla así: 
a o 0) (0) (0) (0) y * 
1. Evaluar g en una aproximación inicial xW) = (a, ii ) . 


2. Determinar una dirección d desde x(%) que origine una disminución del valor de g. 
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3. Desplazar una cantidad adecuada en la dirección d y llamar al nuevo vector xU). 


4. Repetir los pasos del 1 al 3 con el vector x(”. 


Para determinar la dirección de cada x*) desde x"=1 de manera que se garantice una dis- 
minución del valor de g, se introduce la siguiente definición: 


e Definición 6: Si y : R” — R, el gradiente de y en x con x = (1,,12,..., Ty)" se denota 
como Vg(x) y se define: 


vto = (Leo Le... ¿Leo) 


Así la dirección viene dada por un vector s() = —Vg (xl ca notar que se considera el gradi- 
ente negativo dado que se necesita un descenso en los valores de q. 


Luego para determinar el tamaño del desplazamiento, se considera una variable a y una 
función h real-evaluada que cumpla: 


h (a) = 9 (x" + a4s(") 
h! (az) =0 


, donde se toma el valor de a, que produzca el menor valor en h. 


Para obtener el valor de a, primero se suponen a = y a = 1, luego se calculan: 


a = 9 (xa Pvg (0) 


9d = 9 (500 - avg (20) 
(3) 


: (3) (2) (3) 
Si g3 < 91, se acepta a, y se define a)” = =a% O caso O se toma 0, 
se vuelve a realizar el mismo proceso para die un a. Al determinar los tres: 0, 


= ar 1/2 y 
M 


a? y alo y sus respectivos g1 , 92 =y (x*) — al Vg (x0) )) y 93, se considera el EG 


interpolante mediante diferencias divididas de Newton sobre estos tres puntos: 


Pla) = gi +h10a+hgza (a — 0) 


: donde: hi = ELE 
QU — Q1 
has 93 — 92 
QU3g — QU 
ha=h 
hg= =— 
Qg— Q1 


Notar que la solución de P'(a,) = 0 viene dada por: 


Por último se comparan go. = Y (x0) —avy (xM)) con 93, y de nuevo se toma el valor entre 
o. y 03 que genere el menor valor en g y se continúa el proceso. 
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4.4.1 Algoritmo 


x(0) e R” 
x (0) = x (1) + as 4D, k > 1 


; donde s(*=1) = —Vg (x*=D) y az, cumple: 


hilos) =9 (xD) + o -184=0D) 
Aowaj=0 


, y además se toma el a -1 que produzca el menor valor en h. 
4.4.2 Ejemplo 


Aproximar la solución del siguiente sistema con tolerancia de 0.05 y aproximación inicial 
=1(0:0.0+ 
x y) y 


11 +cos (1,1213) — 1 =0 
(1-21) +22 +0.0527—0.1523-1 =0 
—17 — 0.123 +0.01x3 + 23 — 1 ==) 
o Solución: 
Ax) = 2,+c08 (1,1223) — 1 
fx) = (1-21) +22 +0.0522 — 0.1573 — 1 
fax) = —12í-—0.123+0.01%3 + 23 — 1 


(16) + (26) + (869) 
= (2, +c0s (2,22x3) — 1)? + (a —- 21)! + 22 + 0.0523 — 0.1513 — 1) 


5] 
AS 
e 
— 
ll 


+ (—a? 0.107 +0.01x9 + 03 — 1)” 


Vo) = (207 E0O + 200709 +2509 3260, 
2) ZOO +22) 7200 + 210097200) 
ALS + LES +2) E )) 


> Primera iteración: 

g(x0)=1  Vg(x)= (0.1, 0.01, -2) > so = (—0.1,0.01, 2)' 
0% =0.5007 > xY = x0 — aso = (—0.050068074, 0.005006807, 1.001361489)' 
> Segunda iteración: 


g (x"?) =0.0093593  Vg(x”) = (—0.059616, —0.16555, 0.0056453)' 
> s¡ = (0.059616, 0.16555, —0.0056453)' 


01 = 0.5862 > x% = xU — 015, = (—0.050068074, 0.005006807, 1.001361489)' 


> Iteraciones siguientes: 
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Tabla 4: Iteraciones de método de descenso más rápido 


3 xó=1) m0) [E 
0.0000000 —0.050068074 
1 0.0000000 0.005006807 1.001361489 
0.0000000 1.001361489 
—0.050068074 —0.015119550 
2 0.005006807 0.102054203 0.097047396 
1.001361489 0.998052052 
—0.015119550 —0.001526812 
3 0.102054203 0.097176021 0.013592738 
0.998052052 1.000086597 


Se considera la aproximación: x4) = 


ancia de 0.05. 


4.4.3 Código 


(—0.001526812, 0.097176021,1.000086597)* con toler- 


Código 5: Implementación del método de descenso más rápido 


function [x i]= DescensoMR_SENL(Fsis, x0, TOL, n=35) 


AINPUT: 


/Fsis: función de R” a R” asociada al sistema 


40: aprorimación inicial 


ATOL: tolerancia permitida 


/n: numero maximo de iteraciones (35 por defecto) 


AOUTPUT: 


43: aprozimación final 


£i: índice de la ultima ¿iteración 


fprintf('Mtt>> Método de Descenso más rápidon'); 


fprintf('kVtx” (x-1)VtWtx" (x) WWtl | x7 (%)-x" (x-1) | looAn'); 


fprintk(W 


soho LS Ea 


gl = gval(Fsis, x0);%valor de z0 en la función escalar y 


z = ggrad(Fsis, x0, TOL);%valor del gradiende de g en 0 


al = 0.0;%primeros valores de O 


A OA 


v= x0 - a3.*z;%se calcula la prime aprozimación 


g3 = gval(Fsis, v);%se calcula el valor de g en esa aprorimación 


while g3 >= gl/ímientras g3 supere a gl, se buscan nuevos al 


a3 = a3/2; 


v=x0- a3.x*z; 


g3 = gval(Fsis, v); 


endwhile 


a2 = a3/2;/se tiene al y a3, ahora se define al 
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30 


31 


32 


33 


34 


35 


36 


37 


38 


39 


40 


41 


42 


43 


44 


45 


46 


47 


48 


49 


50 


51 


52 


53 


54 


55 


56 


57 


58 


59 


60 


61 


62 


63 


64 


65 


66 


67 


68 


69 


v=xXx0- a2.*z; 


g2 = gval(Fsis, v); 


esto es por las diferencias divididas de Newton 


hi = (g2-81)/a2; 
h2 = (g3-82)/(a3-a2); 
h3 = (h2-h1)/a3; 


ahora se estima el Q 
a0 = 0.5*x(a2-(h1/h3)); 
v=xXx0- a0.*z; 


g0 = gval(Fsis, v); 


if g0 > g3 
a = a3; 
else 
a=a0; 
endif 


x= x0 - a.*z;%se calcula la aprorimación 


err = max(abs(x-x0));%para comparar con la tolerancia (es la norma infinito) 


print ant LO R9EN E MONO ENE LO EN Oda) eri) 


for j= 2:length(x) 


fprintt UN ORINOCO O) DE 


endfor 
fprintt (Y 


if err < TOL %se comprueba sí ya se satisface la tolerancia 


return 


else 


x0 = x; / sí no se satisface, se realiza el mismo proceso 


A/refrescando el valor de 20 


endif 


endfor 


/se llega al marimo de iteraciones 


disp("Numero maximo de iteraciones alcanzado") 


return 


endfunction 


function gx = gval(Fsis,x) 
AINPUT: 

/Fsis: función de IR” a R” 
4x3: vector 

AOUTPUT: 


Código 6: Aproximación de g en x 


asociada al sistema 
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10 


11 


/4gz: valor de gy evaluado en zx; es la suma de los cuadrados 
% de las componentes de F al evaluarse en z. 

total = Fsis(x); 

total = total.”2; 

gx = sum(total); 


endfunction 


Código 7: Aproximación del gradiente de g en x 


function ggx = ggrad(Fsis,x,h) 

AINPUT: 

/Fsis: función de ¡R” a R” asociada al sistema 
Ax: vector 

4h: tamaño de paso 

AOUTPUT: 

/4ggz: vector del gradiente de g en z 


gex = 1:length(x); 


y =X; 

temi = []; 

tem2 = []; 

tem3 = []; 

tem4 = []; 

for i = 1:length(x) 
y(i) = yG)+h; 
teml = Fsis(y); 
tem2 = Fsis(x); 
tem3 = tem1.-tem2; 
tem4 = tem3./h; 
gex(i) = 2*sum(tem4.+tem2); 
y =X,; 

endfor 

ggx = g8x'; 

endfunction 


4.5 Método de continuación 


También se llama método de homotopía e introduce el problema de aproximación dentro de 
una familia de problemas. Específicamente para un problema de tipo F(x) = 0 con solución 
desconocida x*, se considera una familia de problemas descritos mediante un parámetro 
A € [0,1], donde un problema con solución conocida x(0) corresponde a A = 0 y la solución 
desconocida x(1) = x* corresponde a A = 1. 


Bajo estas condiciones se puede definir una función G de esta forma por ejemplo: 
G:[0,1]xR” — R” 
(Ax) — G(A,x) = F(x) + (A — 1)F(x(0)) 


La función G con el parámetro A proporciona una familia de funciones que pueden conducir 
del valor conocido x(0) a la solución x(1) =x*. La función G se llama homotopía entre la 
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función G(0, x) = F(x) —F (x(0)) y la función G(1,x) = F(x). El problema de continuación 
consiste en: Determinar una forma de pasar de la solución conocida x(0) de G(0,x) a la 
solución desconocida x(1) =x* de G(1,x) = F(x). 


El conjunto (x(A)|0 <A < 1) puede verse como una curva parametrizada por Á que va 
desde x(0) =x%) a x(1) = x”. Si las funciones G(A,x) y x(A) son diferenciables, al derivar 
G(A, x) = 0 respecto a Á se tiene que: 


IA KA) IA KA) 
DA | A AS 


Despejando x'(A): 


Es A LOA) 


Como G(A, x) = F(x) + (A — 1)F(x(0)), se tiene que: 


E J(x(A)) y 600) = F(x(0)) , por lo que: 


x (A) = — Ex) F(x(0)) , para 0 <A < 1 y x(0) conocido. 


El siguiente teorema da condiciones bajo las cuales es factible el método de continuación: 


e Teorema 7: Sea F(x) una función continuamente diferenciable en |¡R” para x € KR”. 
Suponga que la matriz Jacobiana Jp(x) es no singular para cada x € R” y que existe una 
constante M tal que || Jp(x)7*] < M, para cada x € R”. Entonces, para cualquier x(0) en 
R”, existe una única solución x(A), tal que: 


G(A,x(A)) =0 
, para cada A € [0,1]. Además, x(A) es continuamente diferenciable y 


x'(A) = — [Jp(x(4))*F(x(0)) , para cada A € [0, 1] 


El sistema de ecuaciones que se debe resolver tiene esta forma: 


dx 
e = OLA Too En) 
dx 
229 q P2lA, T1, La, ..., En) 
dx, 


AA = Dal Bi la) 


onde 
asia A(0)) 
A A O 
dad £u(xc(0)) 
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Para resolverlo se utiliza el método de Runge-Kutta de cuarto orden, primero se elige un 


entero de N >0yh= E = 5. Se divide el intervalo [0,1] en N subintervalos con los 
puntos de red: Aj = ¿h, para cada ¿ = 0,1,2,...,N. Se usa la notación w;; para cada 
¿=0,1,2,...,N y j= 0,1,2,...,n, para denotar una aproximación a 2;(A;). Para la 
condición inicial, se toma: W1 y = 11(0),Wwa0 = 22(0),...,Wno = Tn[0). 
Suponiendo que se han calculado Ww: ;,Ww2;,...,Wn,¡, Se Obtiene W1 ¿+1, W2¿+1) +++, Un ¿+1 US- 
ando las ecuaciones: 
kr; = ho; (As, Wi,j) Wa jj)... ¿Wiy) , para cada 1 = 1, 2, OÍ 
h k k rs 
kai = ho; (a, So 1 72] 7 A ma) ¿para cidad 2 53d 
h k k Kan 
K3 == ho; (a 9 y Wi) > , Wa)j E Wi) j Ea) , para cada ¿1 = 1, 2, ...,N 


Kai = ho; (Aj + h, Wi; + K3,1, Wa; + 32, o , Wn,; + K3.n) , para cada 1= 1, 2, e TO 


y finalmente 


1 
Wij+1 = Wij¡ 6 (5% 2ko 2k3 5 Kka:) , para cada 1 = 1, 2, ...,N 


Se usa la siguiente notación para simplificar la presentación: 


ki1 ka 1 k3,1 kar Wi,j 

k12 ka 2 k3,2 ka 2 Wa j 
k¡ = , , ko= ; , k3= . , k¿= , y Wj¡= 

Kia kan Esa Kan Wn,j 


Notar que x(0) =x(Apo) = Wo, y para cada j =0,1,2,...,n se tiene: 


HA, Wi,;j, Wa,j, e ¿Wi ¿) 
PA Wi Wo cc tDag) E 

k = h AR 2 | =h|=J (159,5. 5)] F(x(0)) 
Pal Aj, W15, W2,j >>>, Wn,j) 


= hl-Jp (w,)]" F(x(0)) 


k 
k> = hide (+3) 
k = hide (+5) 


ky = h[-Jp(w;+k3)] *F(x(0)) 


1 


Finalmente x (A,) = x(1) es la aproximación de x*. 
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4.5.1 Algoritmo 


x(0) e R” 
lA) x (Aj) y 3 (ky y | 2ko PE 2k3 k;) 


; donde dado N > 0, se hace h = 1/N y luego para cada A; se tiene 0 < A¡<1 y A; = jh 
para cada ¿ =0,1,2,..., N, y además: 


k = xy)" 


+ = [a(2opeo) 


2 
.- [coro 
ki = [Jr (x(Ay) +k)]] 7 b 


; con b = —hF(x(0)). 
4.5.2 Ejemplo 
Encontrar la solución del siguiente sistema de ecuaciones con N = 2: 


di -— 2x0 +303+8 =0 
31113 + 2x1 — 512 +8 =0 


o Solución: 


N 2>h 172 0.5 > Ao 0,A1 dog =1 
Sea x(0) = (0,0) > b = —hAF(x(0)) = —0.5(8, 8)' = (-4, -4)' 


> Primera iteración: 


20 0 z 0.05 0 
sl O | loa 
a (0.2186 1. (0.2197 1 — (0.2429 
di a de $ 0.8996 )” A 0.9547 
> x (0.5) = + E (k, + 2k) + 2k3 + k1) = (0.2199375, 0.9060504)' 


> Segunda iteración: 


—18.2405 0.4530 
Jes(05)= | 2.4105 Ea 


ds 0.2430 DE 060749 E 0.2786 e — (0.3321 
de 0.9552 )? 1.0626 )” % e 1.3264 


=> x (1) = x (0.5) + 1 (k, + 2k3 + 2k3 + k4) = (0.5000510, 2.0004611)' 


[Tp (x (0.5) En | —0.055515 —0.0052388 | 


—0.027875  —0.21093 


> Iteraciones: 


Tabla 5: Iteraciones de método de continuación 
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j x (Ay) Xx (Aj41) Belda) sep) lL 
0.0000000 0.2199375 

0 0.0000000 ) 0.9060504 0.9060504 
0.2199375 0.5000510 

1 0.9060504 ) 2.0004611 1.0944107 


4.5.3 Código 


Código 8: Implementación de método de continuación 


function [x i]= Continuacion_SENL(Fsis, 
ZINPUT: 

pFsis: función de IR” a R” asociada al 
4x0: aprorimación inicial 

/N: número de iteraciones 

AOUTPUT: 

43: aprozimación final 


£i: índice de la ultima ¿iteración 


x0, N) 


sistema 


fprintf('Mtt>> Método de ContinuaciónAn'); 


Epia ta inte) INS SAD AA O ONDNOS 


fprint£(' 
h = 1/N; 
b = (Fsis(x0)).*(-h); 
LO ON 


> 
1 


AT 


inví(A);%se calcula la inversa 


wr 
e 
il 


ATx*b; 


A = jacob(Fsis, x0+k1.*(1/2), 107 (-10)); 


> 
RH 
l 


inví(A);%se calcula la inversa 
AIxb; 


Y 
N 
l 


A = jacob(Fsis, x0+k2.*(1/2), 107 (-10)); 


> 
RH 
l 


inví(A);%se calcula la inversa 
AIxb; 


wr 
(0%) 
tl 


A = jacob(Fsis, x0+k3, 107” (-10)); 
inví(A);%se calcula la inversa 
AIxb; 


on > 
5 OH 
Il ñ 


X= Xx0+ (k1 + k2.x2 + k3.x2 + k4) 


err = max(abs(x-x0));%para comparar con la tolerancia (es la norma infinito) 


(1/6); 
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jacob(Fsis, x0, 10”(-10));%se calcula la matriz jacobiama aplicada a TO 


38 


39 


40 


41 


42 


43 


44 


45 


46 


47 


48 


49 


50 


51 


/se refresca la tabla 
print tan O MENO EN cOn a 0 Na ere 
for j = 2:length(x) 

fprintt NON OEA O) 


endfor 
print 


if i == 
return 
else 
240) 528 
endif 
endfor 


endfunction 


25 


3) 


5.1 


9.3 


Problema de aplicación 

Descripción 
Una empresa de teléfonos desea introducir un nuevo producto para competir con otros. 
La empresa planea invertir un millón de dólares para desarrollar este producto. 


El éxito del producto dependerá de la inversión en campaña de marketing y el precio 
final del producto. 


Esta es la fórmula utilizada por el departamento de marketing para estimar las ventas 
del nuevo producto durante el próximo año: 


20000 + 5y/a — 60p 


: donde a: cantidad para invertir en la campaña de marketing 
p: precio del teléfono 


El coste de producción del teléfono es de 100 dólares por unidad 


¿Cómo podría la empresa maximizar sus ganancias para el próximo año? 


Planteamiento 


Beneficios de las ventas: 
(20000 + 5y/a — 60p) p 


Costos totales de producción: 


(20000 + 5y/a — 60p) 100 


Costos de desarrollo: 
1000000 


Costos de marketing: 


Beneficio total: 


S(a, p) = (20000 + 5ya — 60p) (p — 100) — 1000000 — a 


Solución 
Se considera el gráfico de la superficie generada por S(a, p) en la figura 1. 
Se observa que es posible calcular un máximo local. 


Luego se considera el gradiente de $: 


5(p — 100) 737 — 1 ) 


S Ls 
dia ( 20000 + 5/a — 120p + 6000 
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Figura 1: 
S(a, p) = (20000 + 5/a — 60p) (p — 100) — 1000000 — a 


0 


a 


-200000 
J 

-400000 
= 

-600000 
a 
-800000 
-1e+06 
-1.2e+06 


-1.4e+06 
0 


350 


100000 300 


150000 — < 
200 
200000 150 


e Se alcanza el máximo local cuando VS(a,p) =0 


e Se tiene el siguiente sistema: 


e Para determinar una aproximación a su solución se crea el siguiente código en Oc- 
TAVE/MATLAB, donde se utiliza la implementación del método de Newton antes pre- 
sentado (ver código 2 y código 3) 


e Observando el gráfico, se toma una aproximación inicial: x(% = (100000, 100* y tol- 
erancia de 1079, además el código genera el gráfico de la superficie (figura 1) y se 
determina también el valor de la ganancia máxima. 


5(p — 100) 717 — 1 =0 
20000 + 5/a — 120p + 6000 =0 


Código 9: Problema 


HfProblema 
*I) Para el gráfico de la superficie S: 
Para la malla 


[1:1000:200000] ; 
[151:1:350]; 


= meshgrid(x, y); 


Para la superficie 
200) 


= 1:200 
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AE 
MEE 


Z(1,3) = (20000+5*sqrt (a)-60*p) * (p-100)-1000000-a; 


endfor 


endfor 


/Se genera el gráfico de la superficie 


mesh(X,Y,Z) 


II) Para resolver el sistema 


Primera función componente del gradiente 
si = O(a,p) (5* (p-100)*(1/(2*sqrt(a)))-1); 
/Segunda función componente del gradiente 
s2 = Q(a,p) (20000+5x* (sqrt (a) )-120x*p+6000) ; 


Gradiente de S 


Sgrad = 0(x) [s1(x(0),x(2));82(x(1),x(22)]; 


/Aprozrimación inicial 
x0 = [100000;100] ; 


Aplicando implementación 


[x i] = Newton_SENL(Sgrad,x0,107-5,100) 


*III) Para la ganancia máxima 


Ss = Q(a,p) ((20000+5x*sqrt (a) -60x*p) * (p-100)-1000000-a) ; 


SG(1M,0)) 


Tabla 6: Iteraciones 


calculan los valores de Z para cada punto en la malla 


e Los resultados de la implementación se muestran en la tabla 6. 


k Y) x 0%) ¡ES = xD) I 
0.0) 
100000 49124.4947221 
) ( 100 ) ( 226.4911063 ia 


49124.4947221 
226.4911063 


80332.5998434 
228.8351279 


31208.0651213 


80332.5998434 
228.8351279 


101069.2416232 
230.0005069 


20736.6817798 


230.0005069 


105826.0755469 
230.2248042 


| 


4756.8339238 


105826.0755469 
230.2248042 


| 101069.2416232 


) 
) 


( 106003.0165429 


230.2325481724832 


) 


176.9409959 


106003.2449961 
230.2323981 


) 


230.2325981 


106003.0165429 106003.2449961 
' ( 230.2325481724832 ) ( 230.2325581 ) rias 
7 ( 106003.2449973 ) Oia 
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5.4 


La implementación da como resultado una aproximación x(") = (106003.2449973, 230.2325581)' 
con tolerancia de 107? un valor de $ — 88372 como ganancia máxima. 


Análisis de resultados 


Para maximizar las ganancias se deben invertir $106003.24 en marketing y establecer 
un precio de venta de $230.23 por cada teléfono. 


Al sustituir las variables obtenidas en la ecuación S, se estima la ganancia máxima en 


$ — 88372. 


Se concluye que la venta del teléfono no dejará ganancias, sino lo contrario, dejará 
pérdidas y se recomienda no producir el teléfono con esas condiciones. 
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6.1 


6.2 


Conclusiones 


Sobre los métodos de solución 


. El método de punto fijo requiere de varias condiciones que garanticen la convergencia 


de la solución como se ve en el teorema 3, además de que es necesario adaptar el sistema 
para cambiar el problema de búsqueda de solución a un problema de punto fijo como se 
manifestó en el ejemplo correspondiente, sin embargo hay que destacar que no requiere 
de muchos cálculos por cada iteración, el algoritmo correspondiente es bastante simple. 
Pero también hay que mencionar que su velocidad de convergencia es lenta. 


El método de Newton es bastante sólido en velocidad de convergencia, pero requiere de 
condiciones fuertes para garantizar la convergencia, como se vio en el algoritmo corre- 
spondiente a este método es necesario que la matriz Jacobiana de la función asociada 
al sistema sea invertible al aplicarse a cada aproximación, además se realizan bastantes 
Operaciones por cada iteración. 


El método de Broyden relaja la condición de la matriz Jacobiana del método de Newton 
al proponer una aproximación para esas matrices, lo que reduce los cálculos pero pierde 
algo de la velocidad de convergencia que tenía el método de Newton. 


Lo métodos de Newton y Broyden requieren de una buena aproximación inicial dado 
que se basan en un análisis de la matriz Jacobiana que involucra derivadas, se podría 
decir que son "métodos locales”, esto es válido también para los métodos cuasi- Newton 
en general. 


El método de descenso más rápido puede complementar a los métodos cuasi- Newton 
al converger a un punto que sería cercano a la solución verdadera, ese punto de conver- 
gencia podría tomarse como aproximación inicial para los métodos cuasi-Newton. 


El método de descenso más rápido se basa un proceso generalizado para búsqueda de 
puntos críticos , haciendo unas modificaciones en el algoritmo correspondiente como 
por ejemplo tomar el gradiente de g con signo positivo y buscar puntos que maximicen 
a la función g, y así el proceso antes expuesto podría usarse para calcular máximos de 
una superficie. 


En el método de continuación se realizan 4 inversiones por iteración correspondientes al 
cálculo de los vectores k,, k>, kz y k,, en comparación el método de Newton realiza una 
inversión por iteración, una alternativa a esto sería utilizar el método de Runge-Kutta 
de orden 2. 


La implementación de los métodos en el lenguaje de OCTAVE/MATLAB fue relativa- 
mente fácil, notar que para los cálculos que involucran derivadas se ha utilizado un 
aproximación sencilla y no se han utilizado paquetes externos. 


Sobre el problema de aplicación 


. Se utilizó la observación del gráfico correspondiente a la superficie de la ecuación S 


(figura 1) para calcular la aproximación inicial, de no tener esta información se pudo 
optar por utilizar otro método auxiliar para encontrar la aproximación inicial. 
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2. Las observaciones del gráfico fueron importantes, de no tener esa ventaja como en 
sistemas de más de dos variables, al aplicar algún método de solución sobre el gradiente 
de S, no se tendría claridad sobre si el un máximo o un mínimo, lo que implica otro tipo 
de análisis que involucra la segunda derivada, en concreto se puede tomar en cuenta 
la matriz Hessiana de S, concepto que no se ha definido en este material pero que 
determinaría dado un punto crítico si este es máximo o mínimo. 


3. El problema de maximización se transformó en problema de búsqueda de solución 
de un sistema de ecuaciones no lineales, esto también puede usarse para problemas 
de minimización o usar directamente el proceso expuesto el método de descenso más 
rápido. 
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Enlaces 


Método de Broyden: 


— https://ocw.ehu.eus/file.php/81/capitulo-23.pdf 
— https://sites.google.com/site/cm4312014/sis/metodo-de-broyden 


Método de descenso más rápido: 
— https://sites.google.com/site/cm4312014/sis/metodo-dese 
Método de continuación: 


— http://personales.upv.es/dginesta/docencia/posgrado/prac_no_lineal.pdf 


— https: //www.dropbox.com/sc1/fi/aukvccnixw1d3t1g2wvp6/Soluci/C3/B3n-—num?, 
C3/A9rica-de-sistemas-de-ecuaciones-no-lineales.paper?dl=0%r1lkey=4m4kkctov8dkl 


Problema de aplicación: 


— https://rua.ua.es/dspace/bitstream/10045/16373/3/Microsoft/20Word/20-% 
203. SISTEMAS/20ECUACIONES/20N0%20LINEALES . pdf 
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